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Abstract 



We propose Fourier transform (hereafter FT) method for processing images of ex- 
tensive air showers (EAS) detected by imaging atmospheric Cherenkov telescopes 
(IACT) used in the very high energy (VHE) gamma-ray astronomy. The method 
is based on the discrete Fourier transforms (DFT) on compact Lie groups, and the 
use of continuous extension of the inverse discrete transform for approximation of 
discrete EAS images by continuous EAS brightness distribution functions. Here we 
Qh' present the FT-method for SU(3) group. It allows practical realization of the DFT 

method for functions sampled on hexagonal symmetry grids implemented in most 
of the current IACT cameras. We note that the proposed FT-method can also be 
implemented for the rectangular grids using the DFT on SU(2) x SU(2) group. 

We present the results of application of the FT-method to Monte-Carlo simu- 
lated bank of TeV proton and gamma-ray EAS images for a stand-alone telescope. 
Comparing between the FT-method and the currently used standard method for 
signal enhancement, on the basis of parameter ALPHA, shows that the FT tech- 
nique allows a better and systematic increase of the gamma-ray signal. The relative 
difference between these two methods becomes more profound especially for 'pho- 
ton poor' images, for which the standard method deteriorates. It suggests that the 
EAS detection thresholds could be effectively reduced with implementation of the 
FT technique for IACTs. This prediction is further supported by a significant noise 
suppression capability of the method using simple 'low-pass' filters in the image fre- 
quency domain. This new approach allows very deep 'tail' (and 'height') image cuts, 
differentiation of images, frequency spectra, etc., which can be used for development 
of new effective parameters for the EAS image processing. 
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processing 
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1 Introduction 



Since last 15 years, starting from the first confident detection of the Crab Neb- 
ula [1] , the IACTs have proved powerful detectors of very high energy ( VHE, 
conventionally E > 100 GeV) gamma-rays sources. Due to IACTs, and in par- 
ticular, stereoscopic IACT systems (see [2]), VHE gamma-ray astronomy is 
now an established and rapidly growing observational astronomy, with proven 
experimental methodology and with already significant number of detected 
sources of various types (see [3] for the comprehensive recent review). As well 
known, the power of IACTs is based on the realization of the idea [4] that 
extensive air showers (EAS) induced by cosmic rays (CRs) and gamma-rays 
can be effectively discriminated from each other exploiting the intrinsic dif- 
ferences between their Cherenkov light images detected by the multichannel 
optical photoreceivers (cameras) of the telescopes. 

The standard method (hereafter, S-method) for such image discrimination 
involves various "shape" and "orientation" parameters proposed first by Hillas 
[5], which are based on the first and second-order moments of the EAS image 
{Gj}. Recall that the EAS image is a set of numbers Gj of photoelectrons 
registered by the J-th photomultiplier tubes (PMTs) of the IACT camera. 
The integer index J = 1,2, ...,K corresponds to the discrete 2-dimensional 
spatial (angular) coordinates of the PMTs, J = (j, m). In order to improve the 
performance of the 'standard parameterization', different image preprocessing 
techniques have been suggested (e.g. [6,7,8]). Note that all of these image 
filtering methods are based on the image processing directly in the coordinate 
(i.e. image) space. 

Meanwhile, one of the known powerful methods used for signal processing gen- 
erally is the method of discrete Fourier transform (DFT; see e.g. [9,10]). In this 
method the discrete image is transformed into wave number (or frequency, for 
transformation of the time argument) space, in which different 'filters' can be 
applied. Then the filtered image is recovered by the inverse transform. Yet, this 
approach has not been implemented for the EAS images detected by IACTs. 
Perhaps, it could be partly explained by the hexagonal/triangular symmetry 
of the PMT grids implemented in most of the current IACTs, whereas the 
known methods of DFT are developed for images defined on grids of rectan- 
gular symmetry. 

Here we propose a method which allows DFT of 2-dimensional images sampled 
on hexagonal grids. This technique represents a particular case of implemen- 
tation of a general method for Fourier analysis of multidimensional discrete 
functions which uses the so called orbit functions of compact Lie groups as the 
transform basis [11,12,13]. To distinguish it from standard DFT, we abbre- 
viate the DFT on Lie groups as DGT, standing for discrete group transform. 
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Images defined on grids with hexagonal symmetry belong to the case of Fourier 
transform on the group SU(3), whose root lattice displays that symmetry. In 
Section 2 we show that the method is rather simple and can be used in practice 
even without any specific knowledge of the Lie group theory. 

Note that the case of DGT on the group SU(2), which has been considered in 
detail in [14] , is reduced to a specific type of DFT known earlier as discrete 
cosine transform (DCT, see [15]). The 2-dimensional version of DCT is notori- 
ous as the basis for the high-performance image compression standard JEPG. 
It corresponds to DGT on the Lie group SU(2) xSU(2), and allows processing 
of images sampled on the ordinary rectangular grids. Thus, although in this 
paper we describe the Fourier transform method (hereafter, FT-method) for 
the IACTs with hexagonal symmetry of the PMT grids used in most of current 
IACTs [16,17,18,19], the FT-method is equally applicable for IACTs with the 
grids of rectangular symmetry (e.g. [20]). 

Our previous study of DGTs has shown that these transforms are distinguished 
by some valuable properties. In particular, the inverse DGTs can be treated as 
continuous Fourier polynomials, or continuous extensions (CE) of the discrete 
transform which interpolate well the discrete image at any point between the 
grid. Unlike continuous extensions of the standard DFT, but very similar to 
the canonical continuous Fourier transform polynomials, the CEDGTs are 
converging functions. Furthermore, the known properties of localization and 
differentiability of the continuous Fourier transform series hold for CEDGT 
polynomials as well (see [14]). 

After describing in Section 2 the basic mathematical technique of the Fourier 
transform on the Lie group SU(3) and demonstrating the ability of CEDGT 
to approximate the original analytic functions, in Section 3 we apply this tech- 
nique for interpolation of Monte-Carlo (M-C) simulated images of gamma-ray 
and proton induced EAS. The goal of this paper is first of all to present the 
Fourier transform method itself and to discuss some of the related new possi- 
bilities for the IACT image processing, rather than achieving the best result 
for the contemporary stereoscopic IACT arrays. Therefore here we use EAS 
images simulated for a stand-alone IACT with parameters for both the cam- 
era and the collector corresponding to those of the HEGRA telescopes (see 
[22]). It assumes a hexagonal 271-PMT camera with angular size of PMTs 
h = 0.25°, and the telescope collector area ~ 8.5 m 2 . These relatively modest 
parameters result in comparatively coarse ('poor') images. For such images the 
performance of the proposed FT-method, which consists in Fourier transfor- 
mation and then interpolation of discrete images by CEDGT, is significantly 
better than the performance of the S-method. 

In Section 4 we discuss some possibilities for image processing in the fre- 
quency domain that can be opened with the use of the FT-method. These 



3 



possibilities include 'denoising' of the image by cutting off the highest fre- 
quency harmonics before the use of CEDGT. This is a very simple low-pass 
filtering procedure which, nevertheless, appears rather effective for practical 
purposes. We will also demonstrate here, although only qualitatively rather 
than quantitatively (which would require a separate detailed study) that the 
proposed "functional" approach contains a significant potential for reducing 
the effective energy threshold of the IACTs. 

Throughout this paper the comparing with the S-method is done only on 
the basis of parameter ALPHA which has proved to be a very efficient sin- 
gle parameter used in the 'standard parameterization'. In order to keep this 
comparison straightforward, in this paper we refrain from introducing new pa- 
rameters for image discrimination, which in principle would be possible having 
a continuous function instead of a discrete image. 



2 Fourier transforms on SU(3) 

2.1 SU(3) Fourier transform on a grid 

In general terms, the concept 'Fourier transform on a grid' means that a 
discrete function {G k \ k — 0, 1, K} defined on the vector points {r k | k = 
0, 1, K} of a grid in IR 2 (or generally in any n-dimensional real space R"), 
can be represented as a linear combination of wave functions e^ kl ' r - ) sampled 
on the same grid. The set of wave numbers (or 'frequencies' in case of time 
variable) involved is finite, i.e. {kj | i — 0, 1, ...K'}, where K' is not necessarily 
equal to K. It is important, however, that the action of the transform operator 
T would produce an unambiguous 'image' of {G k } in the wave vector space, 
T : {Gk} — > {Aj | j = 0, l,...,Ki} where K\ = K. Since the transform 
operator is linear, the numbers of independent elements in the original and 
transformed images should be equal. 

The method of DFT on compact Lie groups, or the DGT, is based on two 
principal ingredients suggested in [11,12] which allow practical computation 
of Fourier transforms on Lie groups. These are the discretization of the space 
of variables in accordance with the elements of finite order (EFO) of the 
group, and the use of orbit functions of the group as the transform basis. 
The basic property that makes the method work, is the orthogonality of the 
orbit functions on the sets of EFOs (see [11,12,13] for details). The EFOs 
of adjoint orders dividing an integer N make a discrete grid of a definite 
symmetry (depending on the group) in the so called fundamental region of 
the group. 
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Without going into specific mathematical details, which are being presented 
elsewhere, but in order to show the connection between the transform basis and 
the Lie groups, we remind that the group SU(3) can be faithfully represented 
as the group of all unitary matrices of size 3x3. Any unitary matrix can be 
diagonalized by a unitary transformation. Therefore every element of SU(3) 
is conjugate to at least one diagonal element in this 3-dimensional matrix 
representation. The set of all such diagonal matrices forms the maximal torus 
T of the group: 



T = { 



Q(x,y,z) 
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x + y + z = 



The condition x + y + z = describes a plane on which the unitary condition 
det = 1 for the matrices is satisfied. Every element of the group can be thus 
represented by a 2-dimensional vector point on this plane. In particular, the 
simple roots of SU(3), given by vectors «i = (1,-1,0) and a 2 = (0,1,-1), 
correspond to the 2-dimensional vectors with lengths equal to \pl at an angle 
120° to each other. In that plane, the maximal torus T of the group corresponds 
to the hexagonal region shown in Fig. la. 

Two elements in T are conjugate, if and only if they differ by permutation of 
their diagonal entries. The conjugacy classes in SU(3) are in 1-1 correspon- 
dence with the points in the fundamental domain F of the group. In Fig. 1 
it corresponds to the triangular segment enclosed between vectors oj\ and uj 2 
representing the fundamental weights of SU(3). These vectors are expressed 
through aii and a 2 as u\ = (2ai + a 2 )/3 and u 2 = («i + 2a 2 )/3. The weights 
Ui and u 2 are orthogonal to a 2 and a±, respectively, with the scalar product 
{uji ,aj) = 5ij, where 5ij is the Kroneker symbol. 

Any vector A in the plane can be expressed through fundamental weights as 
A = auoi + buo 2 . Action of the Weyl group W on A, which in case of SU(3) 
group is reduced to reflections of A in the lines of vectors lo\ and u 2 , produces 
a finite set {VE(A)} of vectors equidistant from the origin. This set represents 
the Weyl group orbit of A containing the following 6 elements /i: 

W(\) = {auj\ + buj 2 , bui — (a + b)uj 2 , —(a + b)u>i + auj 2 , 

— auj\ + (a + b)u 2 , (a + b)u>i — buj 2 , —buj^ — au 2 } (2) 

In particular cases, when either a = or b = 0, the number of different 
elements in Eq.(2) is reduced from 6 to 3 (see Fig. la), and it is only 1 if 
a = b = 0. Note that the action of the Weyl group expands the fundamental 
region F onto the maximal torus T. 



5 



An orbit function ty/^) = ^a,b at the vector point r (hereafter we use this 
simple notation instead of r) is defined [11,12] as a finite sum of exponential 
functions on the Weyl group orbit of the element A: 

<Mr)= E (3) 



where (//, r) is the scalar product of 2 vectors. Here we will be interested in 
orbit functions corresponding to integer indices a and b. 

An EFO is an element X of the group, such that X N is the unit element E 
for some natural number N. In the cu-basis representation X = a\LOi + b\U02 
this condition is satisfied if (and only if) a\ and b\ are rational numbers. 
The sets of all EFOs of adjoint order dividing N are represented by vector 
points X — > r k ^ m with a\ = k/N and b\ = m/N that satisfy the conditions 
< k,m,k + m < N. These make sets {r km } = F N C F in the form of 
equilateral triangular grids, as shown in Fig. lb for N = 12. 

The key property for the method of DFT on compact semisimple Lie groups 
is the discrete orthogonality of orbit functions with different A on the sets 
Fn of EFOs of adjoint order N [11,12]. In case of SU(3) group this property 
corresponds to the equation 

E Pkm^j,n( r km)^i,p{ r km) = 5i,jSp, n D N (j,n) (4) 



where 



km (l + <Wm,o)(l+4,0 + <W)) ' ^ 



is a multiplicity factor 1 , and the Kroneker symbols are treated modulo N, i.e. 
Si j = 1 if % = j mod(A), e.g. 5n,o = 1- The normalization factor in Eq.(4) is 

^ , x 108A 2 

D N (j,n) = ——, (6) 

-r jn 



Important statements for the method of DFT on Lie groups [11,12] are the 
following two: 

(a) for any given N the set of orbit functions \&_,- jn with indices < j,n, j + n < 
N makes a full set of functions orthogonal to each other on the equilateral 

1 it shows the number of elements in the torus T which are conjugate to the given 
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triangular grid {r km | < k, m, k + m < N} in the form of Eq.(4); 

(b) this set provides a basis for DGT on SU(3) orbit functions with the smallest 

possible wave numbers (j, n). 

A DGT of a given discrete function {G km = G(r km )} produced by sampling 
of a continuous (analog) function G(r) at the grid points {rkm} corresponds 
to solution of the set of equations 

j+n<N 

G km = A jn^jn(r km ) for < k, m, k + m < N . (7) 

j,n>0 



Multiplying Eq.(7) by P km ^ P s( r km), then summing up over {k, m}, and using 
Eq.(4), the transform values {Aj n } are readily found: 



^ k+m<N 



Ajn — n 7 . r X! PkmGkm^jn(fkm)- (8) 



2.2 Continuous extension of DGT 



The set {A km \ < k, m, k + m < N} is an exact Fourier transform of {Gk m }, 
as far as the knowledge of {A km } allows unambiguous reconstruction of {G km } 
at the points {r km } of the given grid using the inverse DGT, i.e. Eq.(7). 

The inverse DGT can be extended into a continuous function in the form of 
Fourier polynomial (i.e. a series of cosine and sine functions) as suggested in 
[14]. It is done simply by replacing the discrete spatial argument r km of the 
orbit functions ^j n { r km) i n the inverse DGT series, Eq.(7), by the continuous 
argument r. The resulting function 

j+n<N 

F N (r)= £ A jn * jn (r) (9) 

j,n>0 



is called continuous extension of DGT, or CEDGT. Note that in case of real 
functions, F N (r) is reduced to J2{j, n }[^ eA -jn Re^ r j n (r) — ImA,- n Im\F, n (r) ]. 

The function F N (r) coincides, obviously, with {G km } = G(r km ) at the grid 
points. However, the most valuable property of F N (r) is in the quality of its 
interpolation of the values G km between the grid points. The approximation to 
the original G(r) that it provides is such that even the first and (under certain 
conditions) also the second derivatives of F^(r) are meaningful functions, con- 
verging to the respective derivatives of G(r). The CEDGT also satisfies the 
property of locality of the transform. These are similar to the properties of 
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canonical continuous Fourier transforms. Note that these properties do not 
hold, however, for the continuous extension of the traditional version of DFT 
(i.e. which is often used in case of rectangular grids). Even the straightfor- 
ward continuous extension of the ordinary DFT is a meaningless (profoundly 
oscillating) function (see [14] for details). 

For convenience here we write the explicit expressions for the SU(3) orbit 
function ^/ a 6 (r). In Cartessian coordinates with the y-axis bisecting the equi- 
lateral triangle (i.e. the fundamental region), with the sides of length 1, as in 
Fig. lb, the real and imaginary parts of the orbit functions are reduced to the 
following: 



■n t / \ ( a + b\ ( a 

Re v a;b (x, y) — 2 cos 2ny — ^=- cos 2irx- 



J V 3 / 

/ a \ ( 2b + a\ , . 

+ 2 cos I 2ixy-j= \ cos \2ttx— — J (10) 

/ b \ ( 2a + b\ 
+ 2 cos \ 2ity-j= \ cos \ 2nx— — I , 

, . n / a + b\ . / a — b\ 
Im y a , b (x, y) = 2 cos \27ry—^-J sin \2ttx—^— J 

( a \ ( 2b + a\ , . 

I 2 cos ( 2iiy-j= \ sin I 2-kx — - — J (11) 

/ b \ ( 2a + b\ 

— 2 cos 2-ny—pz sin 2rcx 

V V3j \ 3 J 

In this system of coordinates, the property of complex conjugation of or- 
bit functions ty a> b = ^b,a results in useful symmetry properties of the real 
and imaginary parts with respect to variables x and y. Also, if a = b, then 
Imty ata (x,y) = 0. Note that in case of integer a and b the harmonic order of 
these polynomials is not always integer. 



2.3 Examples of CEDGT for analytic functions sampled on grids 



The interpolation power of continuous extensions of DGTs is demonstrated in 
Fig. 2. For an example, here we interpolate a function sampled on a rectangular 
grid. The appropriate Lie group for this symmetry is SU(2) xSU(2). The DGT 
in this case happens to result in 2-dimensional case of a transform known as 
DCT (discrete cosine transform). The orbit functions of the SU(2) group, 
for a one-dimensional function, are composed only of 2 exponents, ^\(6) = 
e i2ir\e + e -i2n\e = 2cos(27rA#) with 9 E [0, 1/2] [14]. The basis for the two- 
dimensional DCT simply is a product of 2 cosine functions in 2 orthogonal 
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directions in the plane. 

The discrete image shown in Fig. 2a is produced from sampling of the origi- 
nal analytic function G(r) composed of a sum of two 2-dimensional Gaussian 
distributions with large axes inclined at an angle 25° to each other. The char- 
acteristic widths of both are chosen to be smaller than the distance between 
the grid points, and correspond to dispersions <j± = 2/3N. This grid is rather 
sparse for the image, which is apparent on Fig. 2a. The result of application of 
CEDGT is demonstrated in Fig. 2b (panel on the right). The original function 
is recovered almost exactly, clearly separating the 2 ellipsoids and recovering 
their orientations. The only noticeable difference between CEDGT function 
-Fjv(r) and the original G(r) is reduced to some low-amplitude wiggles. In 
Fig. 2b these wiggles result in the appearance of the dashed contour line cor- 
responding to the level of —0.1% of the maximum intensity in the image. 

It should be noted that a straightforward interpolation of discrete images by 
CEDGT can be useful in case of low level of noise in the signal, but it may be 
less effective if the level of noise is high. This is because the CEDGT 'zooming' 
of an image containing random noise makes the latter more apparent at small 
spatial scales. In this regard, an important property is that the Fourier trans- 
form of an additive random ('white') noise also is a random noise, whereas 
the DGT representation {A, n } of any meaningful image is concentrated in the 
domain of long frequency harmonics (e.g. see [15]). That is, the coefficients 
| Aj n } | have a general tendency to decrease with the increase of j and n. 
Therefore the noise typically becomes relatively stronger, and even may com- 
pletely dominate the signal at high frequencies, as demonstrated in Fig. 3. 
The big dots here correspond to the mean absolute values | A^ m | of the trans- 
form coefficients of different harmonic orders K = k + m < N for an image 
presented in Fig. 4a. The stars show the mean values of the DGT coefficients 
of the noise (contributed by 2 'hot pixels' described below in Fig. 4) in this 
image alone. It is obvious that at high frequencies K > N/2 the the DGT 
image is totally dominated by this noise. 

This property can be exploited for organizing a very simple, but also very 
effective, low-pass filter. Namely, we can cut off in the frequency domain all 
modes A^ m of harmonic orders exceeding some large N cut , i.e. putting A^ m — > 
for indices N cut < K = (k + m) < N. Then we can use the truncated CEDGT 
series that contains only the modes K < N cut to approximate the image. This 
filter can be described by a parameter Cf = 1 — N cut /N such that Cf — 
corresponds to "no filtering". 

The result of application of such low-pass filter is shown in Fig. 4. The orig- 
inal image in Fig. 4a corresponds to continuous function in the form of 2- 
dimensional Gaussian distribution with effective width a± = 1/N — 0.05 and 
length cry = 2a±. Onto it 2 isolated 'hot pixels', with heights of each equal to 
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the half of the maximum of the main image, are superimposed. One of the 'hot 
pixels' is chosen far from the main structure, but the second one falls within 
the image close to its maximum. This results in the total signal in that pixel 
exceeding the maximum value of the main function. 

Fig. 4b (panels in the middle) demonstrates the result of straightforward ex- 
tension of the discrete image using CEDGT without filtering. Note that the 
amplitude of oscillations of the hexagonal Fourier waves in the near vicinity 
of the isolated 'hot pixel' may reach up to ~ 15-20 % of the 'hot pixel' value. 
The distance between the contour lines corresponds to 5 % of the maximum 
of the ellipsoid, and dashed contours correspond to the negative values of the 
CEDGT image. Note that the amplitude of the oscillations rapidly declines 
with distance from the 'hot pixel' because of the localization property of the 
CEDGT [14]. 

In Fig. 4c we show the result of application of the low-pass filter with the 
filtering parameter Cf — 0.5. The main signal practically does not change. 
This is because for the original images with effective widths a± > 1/N, the 
Fourier transformed 'image' {Ak m } is concentrated mostly in the domain of 
low wave numbers (k + m) < N/2, as shown in Fig. 3. Note that even for more 
narrow images, with widths about only one half of the size of pixels in the grid, 
<t_i_ ~ 1/2N, the amplitude of the image after filtering with Cf — 0.5 drops 
only by ~ 25-30 %, but the orientation of the original images is recovered 
quite well. Meanwhile, both "noisy pixels" have almost disappeared from the 
filtered image in the Figure 4c. The signal in the isolated pixel has dropped 
by a factor ~ 4. Even more important is that the intensity in the 'hot pixel' 
inside the main body of the image has practically recovered the expected 
value of the signal. The difference between the contours in Figures 4a and 4c 
is mostly at contour levels corresponding to low intensities. This difference can 
be completely eliminated by modest 'tail' cuts, which is the standard image 
pre-processing procedure discarding the pixels with signals below some level, 
or perhaps even more effectively, by the 'height' cuts as we explain below in 
Section 3. 



3 Performance of the FT method 

To compare relative performances of the FT- and S- methods, Monte-Carlo 
data bank for proton and gamma-ray primaries has been simulated using 
Hillas' MOCCA code [23] complemented with the full ray-tracing of the inci- 
dent Cherenkov photons [24]. The IACT parameters, such as photocollector 
radius/size, altitude, etc., have been assumed corresponding to the telescopes 
of the HEGRA system located at 2200 m above sea level (see [22]). In par- 
ticular, the images are formed in the plane of photoreceiver consisting of 271 
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PMTs (or pixels), with the angular size h = 0.25° arranged in the hexagonal 
grid. This grid is relatively coarse compared to the cameras of the contem- 
porary IACT projects, and also the number of PMTs is relatively modest. 
Both these factors do help, however, to compare and reveal more easily the 
differences between the FT- and S- methods. 

For both proton and gamma-ray events, the showers have been generated 
with energies starting from 0.3TeV, which is significantly lower than the de- 
tection threshold of IACTs with the given parameters. The simulations were 
performed for a telescope pointed to the zenith. The gamma-ray EAS with 
power-law index of the primary photon spectrum —2.6 are generated in the 
vertical direction, i.e. parallel to the telescope optical axis. For the CR pro- 
tons, an isotropic distribution within a vertical cone with the opening (half-) 
angle < 2.7°, and with the spectral index —2.67 has been used. The trigger 
condition for EAS detection was chosen as '2nn/217 > 10pe\ i.e. > 10 pho- 
toelectrons in any 2 near neighbors from 217 inner PMTs. This excludes the 
PMTs of the last (hexagonal) ring from the trigger. In our data bank the 
number of proton EAS passing the trigger is N p>tot = 2354, and N 7ttat = 4229 
for gamma-ray events. 



3.1 CEDGT representation of the EAS images 

In order to apply the DGT, we have to embed the hexagonal grid of the 
camera into the equilateral triangular grid, i.e. into the fundamental domain 
of SU(3), as it is shown in Fig. 5. The camera with 271 pixels contains M = 9 
'hexagonal rings' of PMTs around the central PMT. We formally add one more 
ring assuming Skm = for the numbers of photoelectrons there. It gives some 
additional space for processing images that may contain non-zero numbers 
of photoelectrons in the 9-th ring of the camera. All this results in a grid 
Fjy having N = 30 equal intervals, and N + 1 = 31 grid points, along the 2 
principal non-orthogonal axes lo\ and uii of the triangle. 

Our next step is to reconstruct the values Gkm for the photoelectron distribu- 
tion (the Cherenkov light brightness) function G(r) at the grid points r km . This 
is a useful pre-processing procedure, as far as the numbers Sk m in the detected 
image represent the values of photoelectrons integrated over the surface of each 
individual PMT. Therefore Skm represents only the mean of the continuous 
distribution function over the PMT, and only in the 0-th order approximation 
it represents the real brightness distribution G(r- km ) = Gkm — Sk m /s pix in the 
centres of PMTs. Here s pix is the PMT surface in square angular units, which 
can be chosen as a normalization factor, i.e. s P i x = 1. 

The values G km are better recovered if we assume that G(r) is a continuous 
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function that can be approximated within each PMT by 2-dimensional Taylor 
series around its centre r km . Integrating over the hexagonal surface of the 
PMT, we arrive at an expression that connects the set {G km ) with {Sk m } 
through the second-order derivatives of G(r) at the grid point r km : 

G km ~ -tUL [G >> (rkm) + G'U &0 {r km ) + G:_ m (r km )] . (12) 



Here G", G" +60 (rh m ) and G"_ 60 (rfc m ) are the second derivatives of G(r) along 
the x-axis and along the directions of fundamental weights uj\ and uo 2 of the 
equilateral triangular grid, respectively, i.e. at the angles ±60° relative to the 
x-axis (see Fig. lb). This equation allows fast iterations, approximating ini- 
tially the derivatives of the grid functions in the standard difference scheme 
approach as f"(r ) = [f (r + h) + / (r — h) — 2f (r )]/ h 2 . Here the notations +h 
and —h indicate the positive and negative shifts, respectively, in the given di- 
rection along the grid from the grid point r . For the subsequent iterations the 
derivatives are calculated using directly the CEDGT function. This approach 
provides a good agreement between the integrals of the brightness distribu- 
tion function and the total numbers {Skm} i n the PMTs already after the 
second iteration. The agreement better than < 1 % for most pixels, except for 
some (but not all) of those with very low numbers S km - The total number of 
photoelectrons in the image is preserved within < 1-2 % accuracy. 

In Fig. 5 we show the images of a E — 1 TeV gamma-ray EAS incident at a 
distance 130m from the telescope. Fig.5a (panel on the left) shows the raw 
image detected by the camera, and Fig. 5b (panel in the middle) shows the 
continuous distribution function provided by CEDGT. Note that the lowest 
contour shown corresponds to the level of 3% from the maximum value of 
the image. It demonstrates that the ripples of the Fourier waves are very 
significantly eliminated already at this very low cutoff level 2 . In Figure 5c the 
image is reconstructed after application of the low-pass filter with parameter 
Cf = 0.35. The image here is unified into an essentially single-core pattern, 
as expected for the gamma-ray events. The contour levels in Figure 5c are 
at the same absolute values as in Figure 5b. Comparing with the contours in 
Figure 5b shows a drop in the amplitude of the maximum by ~ 40%. This is 
because in this particular image the signal in one of the pixels was much higher 
than in all of its nearest neighbors. The filter accepts that as a (statistical) 
fluctuation and smooths it out. The total number of photoelectrons in this 
image makes in fact ~ 85% of the original Stat- 

Figure 6 shows the results of application of the same filter with Cf = 0.35 
to the image of a proton EAS with energy E = 2 TeV falling at a distance 

2 Note that for our M-C data bank the level "3% from the image maximum" cor- 
responds on average to 1 photoelectron. 
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100 m from the telescope. The maximum of the filtered distribution in Fig- 
ure 6c is about 85 % of the unfiltered CEDGT image, and the total number 
of photoelectrons found after integration of the filtered image is about 90% 
of the original Stat- This implies that the effective width of the filtered image 
becomes somewhat larger than of the original. The most important feature 
is, however, that unlike the gamma-ray image, the proton image shows the 
main EAS core but it is still far from unification into a single-core pattern. 
Two distinct 'islands', presumably connected with the EAS pions, are still 
apparent. 

3.2 Application to 'pure' EAS images 

The efficiency of the FT-method for signal enhancement is compared with 
the S-method using the orientation parameter ALPHA. This parameter is 
chosen since it is known as one of the most effective parameters (along with 
AZWIDTH) currently used for the IACT image processing in the standard 
parameterization scheme. It is also known that ALPHA corresponds to the 
deflection angle of the image major axes from the direction to the central pixel 
of the camera viewed from the image centre of mass (i.e. from the presumed 
source direction, see e.g. [21]). The signal enhancements provided by FT- and 
S- methods at different levels of tail-cuts c and ALPHA-'cuts' (i.e. chosing 
ALPHA < a), are compared in terms of Q-factors, Q = Q(c,a) = r) 1 /y/fj p ~. 



are the fractions of gamma-ray and proton EAS, respectively, which remain 
in the pool after corresponding parameter cuts c and a. The maximum values 
Qmax are calculated by applying first a fixed tail-cut, and then varying param- 
eter a. Note that we consider only Q-factors which preserve at least 50% of 
the initial gamma-ray events iV 7)tot , i.e. r^ 7 > 0.5. 

To calculate the FT-image, we first reconstruct the values {G p }, using Eq.(12), 
for the continuous distribution function of photoelectrons in the centres of 
PMTs that would correspond to the total numbers of photoelectrons {S p } in 
the image. Here the subscript p — 1, 2, ...P stands for the grid coordinates 
of each of P pixels involved in the image. The Fourier transform coefficients 
Aij are calculated, and the coefficients which do not pass the given filter Cf 
are discarded from the CEDGT series of Eq.(9). Using this truncated CEDGT 
function, the values of the brightness distribution | k — 1,2, ...,K} with 



Here 



7] 7 = r? 7 (c,a) 



T] p = r} p (c, a) 




(13) 
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K ^> P can be calculated for a new grid with any higher density of points. For 
our calculations we divide each of the 6 equilateral triangles of the hexagon 
(the PMT surface) into 4 equal sub-triangles, and calculate the values Gk in 
the centers of each of these sub-triangles. Thus, the density of points in the 
new image is increased by a factor 24. This also leads to an increase of the 
total number of points in the image approximately (because of some zeros) by 
the same factor, i.e. K ~ 24 P. After that a "tail cut" (or "height cut", see 
Section 3.3) procedure is applied, i.e. Gk — > if Gk < g C ut 

Note that here we implement a slightly modified than the currently used tail 
cut procedure. Namely, the tail cut level g cu t is not fixed to some value that is 
the same for all images, as in the standard procedure. Instead, for each image 
we use the "percentage cut", where the cut level g cut = c x G max is defined 
as a fraction < c < 1 from the image maximum G max to cut, but where 
this fraction is now fixed. This procedure takes better into account that the 
amplitude of possible wiggles in the CEDGT function would linearly increase 
with G m ax- In order to make comparing with the S-method more unambiguous, 
we also use the same "percentage cut" approach for the S-method. We have 
checked that the maximal Q-factors reached in the S-method in case of both 
of these tail-cut schemes are practically the same. Note also that in order to 
reduce the impact of statistical fluctuations due to limited numbers of events 
in the (q,q + Ac; Oij,aj + Aa)-bins, we first apply a standard technique of 
averaging of the numbers of events in the bins in the (c, a) plane over (3 x 3) 
square window of bins centered at the given (i, j)-bin. Only then the Q-factors 
are calculated. This procedure makes all the results more systematic, and 
therefore more conclusive. 

In Figure 7a we present the results of application of the FT-method to our 
M-C data bank. Here we show the maximal values of the Q-factors reached at 
different tail-cuts when the parameter a is allowed to vary. The curves 1 and 2 
show the maximal Q-factors for the S-method and for the FT-method without 
filtering (Cf = 0), respectively. Both provide practically the same Q ma x ~ 
2.92, although at somewhat different tail-cuts. The curve 3 shows the maximal 
Q-factors when the filter with Cf = 0.45 is applied. The absolute maximum of 
that curve is significantly increased, reaching Q max ~ 3.4 In Fig. 7b we show 
the a-behaviours of the Q-factors for the same 3 cases, but when the tail-cuts 
are fixed at those (different) values that produce the respective maximal signal 
enhancements. Systematic gain in the signal enhancement after image filtering 
is apparent in both Figs. 7a and 7b, which proves that the filter does work. 

Remaining in the framework of image discriminators based on standard pa- 
rameterization, it seems reasonable to propose that the FT-method would be 
relatively more effective than the S-method for relatively poor images, when 
the numbers of the pixels involved become small, P ~ 10 and less. This is be- 
cause for 'rich' images containing large numbers of active pixels, the calculated 
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moments of a discrete image should be less affected by the sparsity of the data 
than in case of a coarser image. A photon-'rich' image is sufficiently 'smooth' 
already, therefore an additional smoothing of the image by CEDGT would 
not result in a significantly better reconstruction of the image orientation. It 
would not be so, however, for the photon-poor images. 

The number of active pixels P in the image is correlated with the total number 
of detected photoelectrons S tot . Therefore in order to check the validity of that 
proposition, we divide the data bank into 2 subsets containing images with 
Stot < 200 pe and S to t > 200 pe. The value 200 is chosen from the consideration 
to have 2 subsets containing statistically significant numbers of both gamma- 
ray and proton events in both "photon-poor" and "photon-rich" data banks. 
These numbers are equal to Ni jP = 1321 and N 1>7 = 2843 for the photon-poor 
set, and N 2jP = 1033 and N 2tl = 1386 for the photon-rich set. 

In Figure 8 we compare the behaviours of maximal Q-factors at different 
tail cuts for these 2 subsets separately. For both methods the maximum sig- 
nal enhancement has been increased in case of photon-rich images, reaching 
Qmax ~ 3.8 for both. At the same time, for photon-poor images the difference 
between the standard and FT approach is increased further. It is important 
that for the photon-poor images the deterioration of the maximum Q-factor in 
the FT-method is relatively small, by only ~ 0.2, so Q m ax — 3.24. Meanwhile 
the maximum Q-factor in the S-method declines down to Q ma x = 2.64. 

These results suggest that the Fourier transform method allows to work with 
photon poor images significantly better than S-method does. The signal could 
be distributed only in few pixels, but the orientation of the gamma-ray EAS 
could be still recovered with sufficient accuracy by the FT-method. Note that 
even in hypothetical case of only 2 pixels involved in the image, the method 
attributes a minimum dispersion to it, which is of the order of one half of 
the pixel size. In practice, one can process images containing only > 3 pixels. 
FT-method allows to perform shape analyses involving any deep cuts even for 
these very sparse images. 



3. 3 Processing of noisy images 

Photon-poor images are produced mostly by the events with energies of in- 
cident particles closer to the detection threshold of the telescope. Thus, the 
results shown in Figures 8a and 8b seem to suggest that the FT-method con- 
tains a significant potential to deal with the low-energy events better than 
the S-method. This implies a potential for substantial reduction of the energy 
thresholds of the detectors. The question by how much exactly the IACT en- 
ergy thresholds could be reduced with the use of the proposed here FT-method 
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requires a separate detailed study out of the scope of this paper. 

Here we only note that in practice the number of low-energy events detected 
strongly depends on the hardware trigger condition used. Relaxing the trigger 
condition in order to reduce the detection threshold in real experiments results 
in a very significant increase of the number of events with low signal-to-noise 
ratio. In this regard, an important question is whether the FT-method would 
also allow to deal with images with enhanced level of noise over the entire 
camera that would be copiously detected in case of softening the hardware 
trigger condition. 

In order to address this question, we simulate random 'white' noise over the 
entire camera with the mean number of photoelectrons per pixel Uq = 2. 
Poisson distribution with this mean n has been used to generate the noise 
n p in each of the 271 individual pixels. This noise has been added to the M- 
C images of EAS, resulting in S' { = Si + n^. The trigger condition has been 
applied before adding the noise, in order to have the same EAS images in the 
'pure' and 'noisy' sets. 

For noisy images it can make sense to use 'height cuts', instead of usual 'tail 
cuts'. Using again the fixed percentage cuts g cut = cG max as in the tail cuts, 
the 'height' cuts correspond to image modification procedure G p = max(G p — 
g cu t,fy- After such 'height' cuts the image momenta should be less sensitive 
to the noise in distant pixels than after tail cuts which do not modify the 
signal amplitudes in the pixels passing the cut. Therefore one could expect a 
faster, i.e. at lower cut values c, noise suppression and recovery of the signal at 
height-cuts than at tail-cuts. Note that one should not expect a full recovery of 
the initial Q-factors, because also the pixels containing the signal are affected 
by the noise. 

The recovery of the Q-factors for noisy images by S- and FT-methods is shown 
in Figures 9a and 9b using height cuts and tail cuts, respectively. For compar- 
ison, in both panels we also show the Q-factors for the total set of pure (i.e. 
initial) images without noise. Comparing in these 2 panels the curves with full 
and open squares for the FT- and S-methods, respectively, we can see that 
both types of cuts result practically in the same best values for Q max m each 
of these methods if the images are 'noise-free'. All other curves are for the 
noisy data from which the mean value no of the simulated noise is uniformly 
substructed from each pixel, i.e. S\ = max(S'- — 2; 0). Note that this is similar 
to substruction of a 'pedestal' of 2pe. Both panels demonstrate that the signal 
enhancement capability of the S-method and the FT-method without filtering 
(Cf = 0, curves 1) has been essentially killed by noise. One needs very strong 
height or tail cuts at the level of 30% in order to bring the Q-factors back at 
least to the level of ~ 2. Meanwhile, after cutting off the high frequencies the 
FT-method is able to recover the Q-factors to the level of 2.7. One can also 
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see from comparing the curves marked '2' on both panels that the height cuts 
do provide a faster recovery of the Q-factors than the tail cuts, although the 
final results at very large cuts are similar (as expected). 

In Figs. 10a and 10b we show the rates of recovery of Q-factors when different 
levels of 'pedestal substruction' are applied prior to applying the height cuts. 
Curves marked by numbers k = 2, 4, 6 correspond to noisy data from which 
k photoelectrons are substructed, i.e. for each pixel we take the maximum 
S\ k) = max(S{ - Jfe;0). For the FT-method we use the filter C f = 0.45. The 
maximal Q-factors reached in each of these 3 cases by the FT-method are 
practically the same. The only difference between these cases consists in faster 
recovery of the signal when larger k is substructed, as one expects. Meanwhile, 
with the standard method the maximum attainable Q-factors are significantly 
lower than with the FT-method, and also are rather sensitive to the amount 
of pedestal substruction applied. 

Fast recovery of the Q-factor using FT method for image processing makes 
possible extraction of the signal at lower cut levels. It allows to keep more 
images in the raw data, in the sense that photon-'poor' images can be pro- 
cessed more effectively, and using height or tail image cuts at much deeper 
levels than the standard method would allow. 



4 Discussion 

The technique of discrete Fourier transform on orbit functions of Lie group 
SU(3) allows processing of discrete data/images given on the uniform grids of 
triangular or hexagonal symmetry. This is the symmetry implemented in all of 
the contemporary IACTs, including HESS [16], MAGIC [17], VERITAS [18] 
and CANGAROO-III [19]. In case of rectangular PMT grids, such as used 
in CANGAROO-II [20], the appropriate DGT for implementation of the FT 
approach proposed in this paper corresponds to the group SU(2) [14]. 

The DGTs allow straightforward continuous extensions, in the form of trigono- 
metric polynomials, from the discrete points of the grid to any point on the 
image plane. The CEDGT functions F N {r) are distinguished by a number of 
valuable properties, including (a) convergence to the original function G(r) 
with increasing N; (b) localization property ; and (c) differentiability of i*V(r), 
implying convergence of the derivative of Fjv(r) to G(r) [14]. These proper- 
ties are very similar to the properties of the canonical Fourier transforms of 
continuous functions. 

The analysis on the M-C simulated data bank presented here shows that 
CEDGT approach can be effective for enhancement of gamma-ray signals. Re- 



17 



maining in the framework of standard parameterization, i.e. without exploring 
new possibilities that can be opened with the use of FT approach, the Q-factor 
for the parameter ALPHA is steadily enhanced from Q max approx2.9 in the S- 
method to Q ma x = 3.4 simply by using the CEDGT with the low-pass filter 
Cf — 0.45. The latter is a simple possibility provided by the Fourier method 
to suppress the statistical fluctuations or noise (of different origins) from the 
data. Our calculations indicate that the useful range for the filter parameter 
is within Cf ~ (0.3 — 0.5), depending on the type of image analysis which is 
being carried out. For the orientational parameter ALPHA values closer to 0.5 
seem to be optimal, and in this paper we have fixed it to Cf = 0.45. However, 
one might also chose values ~ 0.3 or even less for studies where a lesser degree 
of broadening and unification of proton images would be advantageous. 

Compared with the standard method, the power of continuous extension of 
DGTs for image processing becomes more apparent for photon-poor images, as 
demonstrated in Figs. 8a and 8b. Furthermore, the FT-method offers to deal 
much better that S-method with images where the level of random noise over 
camera is substantially enhanced. The implication is that perhaps with the 
use of FT-method the effective energy threshold of EAS detected by IACTs 
could be further decreased. A more definite answer to this question requires, 
obviously, separate studies involving relaxed trigger conditions, and possibly 
also considering lower than currently levels of the pedestal substruction. 

Out of the scope of this paper also remain possibilities of the FT-method for 
the stereoscopic systems of IACTs. The single notice in this respect is that 
the power of FT-method for better reconstruction of the EAS direction and 
keeping more gamma-ray events in the "single telescope pool" should obviously 
show up in the systems of IACTs as well (and may only be amplified). 

Neither we consider here any new possibilities for image processing that might 
be possible with this technique. Those methods can be connected either with 
image processing in the image domain itself (such as image analyses under deep 
'tail cuts' and 'frequency cuts', use of image derivatives, etc), or image analyses 
in the frequency domain (analysis and optimization of the FT coefficients). 

At last, we would like to note that the technique of Fourier transforms de- 
scribed here is unique for grids with hexagonal symmetry, and it is easy to 
use in practical calculations. Therefore it might be useful for data processing 
not only for the IACTs, but also for other detectors using hexagonal symme- 
try, such as HIRES detector in ultra-high energy CR physics, or the IceCube 
detector in VHE neutrino astronomy. 
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Fig. 1. a (left panel): Representation of the maximal torus T (hexagon) of SU(3) 
group in the plane of simple roots ct\ and ai of the group, the fundamental domain 
F and the fundamental weights uj\ and 002- Also shown are the Weyl group orbits of 
two elements in F, which produce Weyl groups consisting of 6 elements (shown by 
stars) or 3 elements (big dots); b (right panel): A set of elements of finite adjoint 
order N = 12 (big dots) in the fundamental domain (the triangle). 




Fig. 2. a (left panel) The discrete image produced from sampling of a continuous 
function G(r) composed of two Gaussian ellipsoids with large axes inclined at an 
angle 25° to each other, and with widths corresponding to dispersions a± = 2/3N 
for a rectangular grid with N = 20; b (right panel): The reconstructed CEDGT 
image. The contours correspond to intensity levels separated by 10% of the peak 
intensities in the ellipsoids, starting from the dashed contour at a small negative 
level of -0.1%. 
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Harmonic order K 

Fig. 3. The averaged absolute values | Ak m I of the DGT coefficients of the given 
harmonic order K = k + m (shown by big dots) calculated for the image presented 
in Fig. 4. The stars show the mean values | A' k | of the transform of the 2 'hot 
pixels' (see the caption in Fig. 4) alone. 
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Fig. 4. a (2 panels on the left): The contour plot (upper panel) and the 3-dimensional 
view of an analog image composed of 2-dimensional continuous Gaussian ellipsoid, 
onto which 2 narrow spikes which produce 2 'hot pixels' (after sampling of the 
analog signal on the grid) are superimposed. The heights of the spikes are equal to 
half of the height of the ellipsoid; b (2 panels in the middle): A direct (i.e. without 
filtering) CEDGT interpolation of the discrete image produced by sampling of the 
analog image on the grid with N = 20; c {2 panels on the right): CEDGT view of 
the image after application of the filter with Cf = 0.5 (see text). 
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Fig. 5. a (left panel): The discrete image of a gamma-ray event with energy 
E = 1 TeV falling at the distance 130 m from the telescope; b The contour plot of 
the same event in the CEDGT representation without filtering, Cf = 0. The lowest 
contour shown corresponds to the level of 3 % from the maximum intensity of the 
image. The camera embedded into the fundamental triangle is also shown; c The 
image of the same EAS after application of the low-pass filter with the parameter 
Cj = 0.35 (see text). 




Fig. 6. The same as in Figure 5, but for the proton EAS with energy 2 TeV incident 
at a distance 100 m from the telescope. 
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Tail Cut (%) ALPHA (deg) 



Fig. 7. a (left panel): The maximal Q-factors attained at different tail-cuts in the 
standard approach (curve 1, open squares), and by the Fourier transform method 
without filtering (curve 2, triangles), and in case of low-pass filter with the parameter 
Cf = 0.45 (curve 3, full squares); b (right panel): The depences of the Q-factors 
on a for the same 3 cases, but when the tail-cuts are fixed at the respective values 
corresponding to the absolute maxima of Q m ax of the curves on the left panel. 
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Fig. 8. The maximal Q-factors obtained by the S-method (open squares) and by 
the FT-method using the filter C/ = 0.45, calculated for 2 subsets of the original 
M-C image bank (used for Figure 7) that contain only photon-poor images with 
Stot < 200 (a - left panel), or only photon-rich images with Stat > 200 (b - right 
panel) . 
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Fig. 9. The recovery of the Q-factors with increasing levels of height cuts (a - left 
panel) or tail cuts (b - right panel) for the noisy data sets in case of the S-method 
(open circles), and the FT-method (full circles) without filtering (curve 1) and 
using the filter with Cf = 0.45 (curve 2). For comparison, the open and full squares 
show the Q-factors for the initial (pure) image set. For the noisy images pedestal 
substruction equal to 2pe has been used (see text). 




Fig. 10. The extents of recovery of the signals from the noisy data sets with different 
levels of pedestal substruction (see text) in case of FT-method (a - left panel) 
and S-method (b - right panel). For comparison, for both methods also the initial 
Q-factors corresponding to the set of pure images are presented. 
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